/* matrix.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/fwd.hpp>
#include <src/miscellaneous/OS_Basic.h>
/*! \file matrix.h
	 \brief Proxies to boost::numeric::ublas matrix definitions.
*/
#ifndef _OSMatrix_
#define _OSMatrix_
/*! \brief Ublas dense vector
*/
#define UBVector boost::numeric::ublas::vector 
/*! \brief unspecialized cast from OS Array's to Ublas Vector.
	use std::copy workaround.
	Thus it should be overloaded for special type for performance issues.
*/
template <class Array>
UBVector<typename Array::value_type> cast_au(const Array &data) {
	UBVector<typename Array::value_type> result(data.size());
	std::copy(data.begin(),data.end(),result.begin());
	return result;
};
/*! \brief unspecialized cast from Ublas Vector to OS Array's.
	use std::copy workaround.
	Thus it should be overloaded for special type for performance issues.
*/
template <class Array>
Array cast_ua(const UBVector<typename Array::value_type> &data) {
	Array result(data.size());
	std::copy(data.begin(),data.end(),result.begin());
	return result;
};
/*! \brief default sparse matrix container
*/
#define OSMatrix_sparse boost::numeric::ublas::mapped_matrix

/*! \brief default weight set container is default sparse matrix container
*/
#define weight_set_default_def OSMatrix_sparse

/*! \brief boolean equality operator over sparse matrix.
*/
template <class T>
bool operator == (const OSMatrix_sparse<T> &left,const OSMatrix_sparse<T> &right){
	if (left.size1() != right.size1()) return false;
	if (left.size2() != right.size2()) return false;
	if (left.nnz_capacity() != right.nnz_capacity()) return false;
	bool result = true;

	typedef OSMatrix_sparse<T>::const_iterator1 const_iterator1;
	typedef OSMatrix_sparse<T>::const_iterator2 const_iterator2;
	const_iterator1 itl,itr;
	const_iterator2 itl_,itr_;
	for (itl = left.begin1(),itr = right.begin1()
		;itl != left.end1()
		;++itl,++itr) {
			if (itr == right.end1()) return false;
			for (itl_ = itl.begin(),itr_ = itr.begin()
				;itl_ != itl.end()
				;++itl_,++itr_) {
					int breakpoint1 = itl_.index1();
					int breakpoint2 = itl_.index2();
					if (itr_ == itr.end() || *itr_ != *itl_) return false;
			};
			
	};
	return true;
};
/*! \brief boolean non equality binary operator over sparse matrix.
*/
template <class T>
bool operator != (const OSMatrix_sparse<T> &left,const OSMatrix_sparse<T> &right){
	return !(left==right);
};


/*! \brief right multiplication of vector by a sparse matrix
*/
template <class Array>
const Array operator*(const Array& array_, const OSMatrix_sparse<typename Array::value_type>& m)
{
	typedef  UBVector<typename Array::value_type> UBV ;
     return cast_ua<Array>(prod(cast_au<Array>(array_),m));
};

/*! \brief left multiplication of vector by a sparse matrix
*/
template <class Array>
const Array operator*(const OSMatrix_sparse<typename Array::value_type>& m, const Array& array_) {
		typedef  UBVector<typename Array::value_type> UBV ;
/*		breakpoint
		std::cout << std::endl << m << std::endl;
		std::cout << std::endl << array_ << std::endl;
*/
		return cast_ua<Array>(prod(cast_au<Array>(array_),m));
};

/*! \brief dense matrix definition as Ublas matrix
*/
#define OSMatrix_dense boost::numeric::ublas::matrix


/*! \brief boolean equality binary operator over dense matrix.
*/
template <class T>
bool operator == (const OSMatrix_dense<T> &left,const OSMatrix_dense<T> &right){
	OS_DYNAMIC_CHECK(false," bool operator == (const OSMatrix_sparse : not implemented");
};
/*! \brief boolean non equality binary operator over dense matrix.
*/
template <class T>
bool operator != (const OSMatrix_dense<T> &left,const OSMatrix_dense<T> &right){
	return !(left==right);
};
/*! \brief right multiplication of vector by a dense matrix
*/
template <class Array>
const Array operator*(const Array& array_, const OSMatrix_dense<typename Array::value_type>& m)
{
	typedef  UBVector<typename Array::value_type> UBV ;
     return cast_ua<Array>(prod(cast_au<Array>(array_),m));
};

/*! \brief left multiplication of vector by a dense matrix
*/
template <class Array>
const Array operator*(const OSMatrix_dense<typename Array::value_type>& m, const Array& array_) {
		typedef  UBVector<typename Array::value_type> UBV ;
		return cast_ua<Array>(prod(cast_au<Array>(array_),m));
};

#endif
